% Sameer D. Meshram
% IIT Delhi
% Surface roughness from CT-derived profile image
% Standards-oriented workflow following ISO 16610-21 Gaussian profile filter

clear; close all; clc;

%% --- User inputs ---
% =========================================================================
% Input parameters explained
% =========================================================================
%
% (1) Image filename
%     Path (or file name in the working directory) of the profile image.
%     The image may be grayscale or RGB; it is converted to grayscale and
%     then normalized to uint8 internally, so the threshold in (3) is
%     always on a 0-255 scale regardless of the source dtype.
%     The image should contain a clearly visible boundary between the
%     bright material region and the dark background.
%
% (2) Pixel size (mm/pixel)
%     Physical size of one pixel in millimeters. This is the calibration
%     that turns pixel indices into real distances, so it directly scales
%     ALL roughness results. An error here produces results that are wrong
%     by exactly the same factor.
%     Typical sources: CT scanner metadata (voxel size) or a calibrated
%     scale bar. Examples:
%         5.0 micrometres -> 0.005
%         3.6 micrometres -> 0.0036
%         1.0 micrometres -> 0.001
%
% (3) Threshold for white boundary extraction (0 - 255)
%     Intensity value that separates material (>= threshold) from
%     background (< threshold) after the image is normalized to uint8.
%     A good starting point is 128; increase it if the thresholded mask
%     swallows background speckle, decrease it if parts of the material
%     come out as holes. Inspecting the generated "Extracted edge" figure
%     is the fastest way to tune this value.
%
% (4) Minimum object size (pixels)
%     Connected regions smaller than this many pixels are discarded
%     (bwareaopen). This suppresses isolated noise specks near the surface
%     that would otherwise be picked up as the topmost or bottommost
%     white pixel in a column.
%     Keep it well below the smallest real surface feature you care
%     about. Typical range: 20 - 200 pixels.
%
% (5) Edge to extract: 'top' or 'bottom'
%     Which side of the binary material region becomes the profile.
%     In image coordinates row 1 is at the TOP of the image.
%         'top'    -> topmost white pixel in each column
%         'bottom' -> bottommost white pixel in each column
%     Choose based on which side of the bright region is the physical
%     surface of interest. If the material sits below the bright line
%     in the image (so its surface is the upper edge), use 'top'.
%
% (6) Form removal polynomial degree (0, 1 or 2)
%     Subtracts a least-squares polynomial trend from the profile BEFORE
%     any filtering, separating nominal shape from roughness + waviness.
%         0 -> subtract mean only (keep tilt and curvature)
%         1 -> remove linear tilt  (most common for slightly tilted samples)
%         2 -> remove tilt + quadratic curvature (curved parts: shafts,
%              cylinders, lens surfaces, etc.)
%     Without form removal, a strong tilt bleeds into the Gaussian filter
%     result at the ends and inflates Ra / Rq.
%
% (7) lambda_s (mm, 0 to disable)
%     Short-wavelength cutoff of the S-filter. Removes very short-
%     wavelength content (noise, single-pixel jaggies) below this
%     wavelength before the C-filter is applied. Per ISO 4287 / 16610-21
%     the roughness profile lives between lambda_s and lambda_c.
%     Typical values (see ISO 3274 / 4288 tables):
%         0.0025 mm ( 2.5 um) - very fine finishes
%         0.0080 mm ( 8.0 um) - standard general-purpose setting
%     Set to 0 to disable the S-filter (useful when the data is already
%     clean or when you want Ra / Rq that only depend on lambda_c).
%
% (8) lambda_c (mm)     <-- the single most important parameter
%     Long-wavelength cutoff of the Gaussian profile filter. Separates
%     waviness (wavelengths > lambda_c) from roughness (< lambda_c).
%     Standard ISO 4288 cutoffs, chosen by expected Ra:
%         0.08 mm  for Ra up to ~0.1 um
%         0.25 mm  for Ra ~0.1 - 2 um
%         0.80 mm  for Ra ~2 - 10 um    (commonly used default)
%         2.50 mm  for Ra ~10 - 50 um
%         8.00 mm  for Ra above that
%     The total profile must be long enough to support the evaluation
%     length, i.e. L_total >= n_eval * lambda_c (plus trim, if enabled).
%
% (9) Number of cutoff lengths for the evaluation length (n_eval)
%     The evaluation length used for parameter calculation is
%         le = n_eval * lambda_c.
%     ISO 4287 / 4288 specifies n_eval = 5 (five sampling lengths). The
%     Rz definition used here (mean of peak-to-valley over five sampling
%     lengths of lambda_c) also requires n_eval = 5. Use a smaller value
%     only if the profile is too short for the standard case.
%
% (10) Trim lambda_c/2 from each end of the R-profile (1 = yes, 0 = no)
%     Gaussian filtering introduces systematic edge bias even with mirror
%     padding. ISO practice is to discard the outer lambda_c/2 of the
%     filtered R-profile before placing the evaluation window and
%     computing parameters. Leave this on (1). Turn it off only if the
%     profile is very short and you need every available sample.
% =========================================================================

defaults = { ...
    'R3_a_ana.png', ... % (1)  Image filename (see comment block above)
    '0.0036', ...       % (2)  Pixel size in mm/pixel  (e.g. 0.0036 = 3.6 um). This value has to be calculated for image being analyzed by user.
    '60', ...           % (3)  Threshold 0-255 on the uint8-normalized image. This value requires manual tuning.
    '50', ...           % (4)  Min object size in pixels for speckle removal
    'top', ...          % (5)  Edge to extract: 'top' or 'bottom'
    '1', ...            % (6)  Form removal degree: 0 = mean, 1 = tilt, 2 = curvature
    '0.0025', ...       % (7)  lambda_s in mm (S-filter); 0 disables it
    '0.8', ...          % (8)  lambda_c in mm (Gaussian C-filter, ISO 16610-21)
    '5', ...            % (9)  Number of cutoff lengths:  le = n_eval * lambda_c
    '1' ...             % (10) Trim lambda_c/2 from each end: 1 = on, 0 = off
    };

prompt = { ...
    'Image filename', ...
    'Pixel size (mm/pixel)', ...
    'Threshold for white boundary extraction (0-255 on uint8-normalized image)', ...
    'Minimum object size (pixels)', ...
    'Edge to extract: top or bottom', ...
    'Form removal polynomial degree (0, 1 or 2)', ...
    'lambda_s (mm, 0 to disable S-filter)', ...
    'lambda_c (mm)', ...
    'Number of cutoff lengths for evaluation length', ...
    'Trim lambda_c/2 from each end of R-profile (1 = yes, 0 = no)' ...
    };

dlg_title = 'Surface Roughness Inputs';
answer = inputdlg(prompt, dlg_title, [1 80], defaults);
if isempty(answer)
    error('User cancelled input dialog.');
end

image_file      = strtrim(answer{1});
pixel_size_mm   = str2double(answer{2});
threshold       = str2double(answer{3});
min_object_size = str2double(answer{4});
edge_choice     = lower(strtrim(answer{5}));
form_degree     = str2double(answer{6});
lambda_s_mm     = str2double(answer{7});
lambda_c_mm     = str2double(answer{8});
n_eval          = str2double(answer{9});
trim_edges      = logical(str2double(answer{10}));

validateattributes(pixel_size_mm,   {'numeric'}, {'scalar','real','positive'});
validateattributes(threshold,       {'numeric'}, {'scalar','real','>=',0,'<=',255});
validateattributes(min_object_size, {'numeric'}, {'scalar','real','integer','>=',0});
validateattributes(form_degree,     {'numeric'}, {'scalar','real','integer','>=',0,'<=',2});
validateattributes(lambda_s_mm,     {'numeric'}, {'scalar','real','>=',0});
validateattributes(lambda_c_mm,     {'numeric'}, {'scalar','real','positive'});
validateattributes(n_eval,          {'numeric'}, {'scalar','real','positive'});

if ~ismember(edge_choice, {'top','bottom'})
    error('Edge must be top or bottom.');
end

%% --- Step 1: Read the image and normalize to uint8 ---
img = imread(image_file);
if size(img, 3) == 3
    gray = rgb2gray(img);
else
    gray = img;
end
gray = im2uint8(gray);   % handle any dtype uniformly against a 0-255 threshold

%% --- Step 2: Threshold and clean ---
bw = gray >= threshold;
bw = bwareaopen(bw, min_object_size);

if ~any(bw(:))
    error('Binary image is empty after thresholding. Lower the threshold or check the image.');
end

%% --- Step 3: Per-column edge extraction (one z per x by construction) ---
% Row 1 is at the top of the image. "top" means smallest row index with a
% white pixel in that column; "bottom" means largest row index.

valid_cols = any(bw, 1);
x_pix_all  = find(valid_cols);

switch edge_choice
    case 'top'
        [~, z_pix_top]  = max(bw, [], 1);     % first 1 scanning from the top
        z_pix_all       = z_pix_top(valid_cols);
    case 'bottom'
        bw_flipped      = flipud(bw);
        [~, z_pix_flip] = max(bw_flipped, [], 1);
        z_pix_bot       = size(bw,1) - z_pix_flip + 1;
        z_pix_all       = z_pix_bot(valid_cols);
end

boundaryCoords_pixels = [x_pix_all(:), z_pix_all(:)];

if size(boundaryCoords_pixels,1) < 5
    error('Too few boundary columns extracted. Check threshold or image quality.');
end

%% --- Step 4: Visualize extracted edge on original image ---
figure;
imshow(img); hold on;
plot(boundaryCoords_pixels(:,1), boundaryCoords_pixels(:,2), 'r-', 'LineWidth', 1.5);
title(sprintf('Extracted %s edge', edge_choice));
hold off;

%% --- Step 5: Export pixel coordinates ---
writematrix(boundaryCoords_pixels, 'BoundaryCoordinates_pixels.xlsx');
disp('Boundary coordinates exported to BoundaryCoordinates_pixels.xlsx');

%% --- Step 6: Convert to mm, flip so peaks point up, resample uniformly ---
x_raw = boundaryCoords_pixels(:,1) * pixel_size_mm;
z_raw = boundaryCoords_pixels(:,2) * pixel_size_mm;

% Flip so that physical peaks point upward (image rows grow downward).
z_raw = -z_raw;

% Sort and dedupe (extraction already gives unique x, this is belt-and-braces).
[x_raw, sortIdx] = sort(x_raw);
z_raw = z_raw(sortIdx);
[x_raw, uniqIdx] = unique(x_raw, 'stable');
z_raw = z_raw(uniqIdx);

% Uniform grid at the pixel spacing.
dx = pixel_size_mm;
x  = (x_raw(1):dx:x_raw(end)).';
if numel(x) < 5
    error('Profile too short for analysis.');
end
z  = interp1(x_raw, z_raw, x, 'linear');

L_total = x(end) - x(1);

%% --- Step 7: Form removal (tilt / curvature) ---
if form_degree == 0
    z_form = z - mean(z);
    form_component = mean(z) * ones(size(z));
else
    p = polyfit(x, z, form_degree);
    form_component = polyval(p, x);
    z_form = z - form_component;
end

%% --- Step 8: Build Gaussian kernels (ISO 16610-21) ---
if lambda_s_mm > 0
    g_s = isoGaussianKernel(lambda_s_mm, dx, numel(z_form));
else
    g_s = [];
end
g_c = isoGaussianKernel(lambda_c_mm, dx, numel(z_form));

%% --- Step 9: Raw-profile roughness (mean-centered, no filtering) ---
z_raw_centered = z - mean(z);
Ra_raw = trapz(x, abs(z_raw_centered)) / L_total;
Rq_raw = sqrt(trapz(x, z_raw_centered.^2) / L_total);

%% --- Step 10: Apply optional S-filter ---
if lambda_s_mm > 0
    z_s = mirrorConv(z_form, g_s);
else
    z_s = z_form;
end
z_s_centered = z_s - mean(z_s);
Ra_s = trapz(x, abs(z_s_centered)) / L_total;
Rq_s = sqrt(trapz(x, z_s_centered.^2) / L_total);

%% --- Step 11: Gaussian C-filter and R-profile ---
long_wave  = mirrorConv(z_s, g_c);
r_profile  = z_s - long_wave;
r_profile  = r_profile - mean(r_profile);

Ra_filt_full = trapz(x, abs(r_profile)) / L_total;
Rq_filt_full = sqrt(trapz(x, r_profile.^2) / L_total);

%% --- Step 12: Trim edge bias (optional) ---
if trim_edges
    trim = lambda_c_mm / 2;
    if (L_total - 2*trim) <= 0
        warning('Profile too short to trim lambda_c/2 on each end. Skipping trim.');
        x_valid = x;
        r_valid = r_profile;
    else
        keep    = (x >= x(1) + trim) & (x <= x(end) - trim);
        x_valid = x(keep);
        r_valid = r_profile(keep);
    end
else
    x_valid = x;
    r_valid = r_profile;
end

L_valid = x_valid(end) - x_valid(1);

%% --- Step 13: ISO-style evaluation length subset (centered in valid region) ---
le_target = n_eval * lambda_c_mm;

if le_target > L_valid
    warning(['Requested evaluation length le = %.4f mm exceeds the usable ', ...
             'profile length %.4f mm after trimming. Using full usable length.'], ...
            le_target, L_valid);
    x_eval = x_valid;
    r_eval = r_valid;
    le_used = L_valid;
else
    x_mid   = 0.5 * (x_valid(1) + x_valid(end));
    x_start = x_mid - le_target/2;
    x_end   = x_mid + le_target/2;
    mask    = (x_valid >= x_start) & (x_valid <= x_end);
    x_eval  = x_valid(mask);
    r_eval  = r_valid(mask);

    if numel(x_eval) < 5
        warning('Too few points in target evaluation length. Using full usable length.');
        x_eval = x_valid;
        r_eval = r_valid;
        le_used = L_valid;
    else
        le_used = x_eval(end) - x_eval(1);
    end
end

%% --- Step 14: Roughness parameters on the evaluation length ---
r_eval = r_eval - mean(r_eval);     % final centering for parameter calculation

Ra_eval = trapz(x_eval, abs(r_eval)) / le_used;
Rq_eval = sqrt(trapz(x_eval, r_eval.^2) / le_used);
Rp      = max(r_eval);
Rv      = -min(r_eval);
Rt      = Rp + Rv;

% Rz (ISO 4287): mean of Rz values over 5 consecutive sampling lengths of length lambda_c.
Rz = NaN;
n_samp = 5;
if le_used >= n_samp * lambda_c_mm * (1 - 1e-9)
    Rz_vals = zeros(n_samp,1);
    x_s0 = x_eval(1);
    for k = 1:n_samp
        a = x_s0 + (k-1) * lambda_c_mm;
        b = x_s0 +  k    * lambda_c_mm;
        seg = r_eval((x_eval >= a) & (x_eval <= b));
        if isempty(seg)
            Rz_vals(k) = NaN;
        else
            Rz_vals(k) = max(seg) - min(seg);
        end
    end
    Rz = mean(Rz_vals, 'omitnan');
end

% Skewness and kurtosis relative to Rq.
if Rq_eval > 0
    Rsk = mean(r_eval.^3) / Rq_eval^3;
    Rku = mean(r_eval.^4) / Rq_eval^4;
else
    Rsk = NaN; Rku = NaN;
end

%% --- Step 15: Console summary ---
fprintf('\n--- Surface Roughness Results ---\n');
fprintf('Image file                 : %s\n', image_file);
fprintf('Pixel size                 : %.6f mm/pixel\n', pixel_size_mm);
fprintf('Edge extracted             : %s\n', edge_choice);
fprintf('Form removal polynomial    : degree %d\n', form_degree);
fprintf('Threshold                  : %g (on uint8-normalized image)\n', threshold);
fprintf('Minimum object size        : %g pixels\n', min_object_size);
fprintf('lambda_s                   : %.6f mm\n', lambda_s_mm);
fprintf('lambda_c                   : %.6f mm\n', lambda_c_mm);
fprintf('Edge trim (lambda_c/2)     : %s\n', ternary(trim_edges,'on','off'));
fprintf('Requested eval length      : %.6f mm\n', le_target);
fprintf('Used eval length           : %.6f mm\n', le_used);
fprintf('Total profile length       : %.6f mm\n\n', L_total);

fprintf('Raw mean-centered profile  : Ra = %.6f mm, Rq = %.6f mm\n', Ra_raw, Rq_raw);
fprintf('After form removal + S     : Ra = %.6f mm, Rq = %.6f mm\n', Ra_s,   Rq_s);
fprintf('R-profile, full length     : Ra = %.6f mm, Rq = %.6f mm\n', Ra_filt_full, Rq_filt_full);
fprintf('R-profile, eval length     : Ra = %.6f mm, Rq = %.6f mm\n', Ra_eval, Rq_eval);
fprintf('Rp (max peak)              : %.6f mm\n', Rp);
fprintf('Rv (max valley)            : %.6f mm\n', Rv);
fprintf('Rt (peak-to-valley)        : %.6f mm\n', Rt);
if ~isnan(Rz)
    fprintf('Rz (ISO 4287, 5 samples)   : %.6f mm\n', Rz);
else
    fprintf('Rz (ISO 4287, 5 samples)   : not computed (le_used < 5*lambda_c)\n');
end
fprintf('Rsk (skewness)             : %.6f\n', Rsk);
fprintf('Rku (kurtosis)             : %.6f\n', Rku);

%% --- Step 16: Plots ---
% Plot 1: profile filtering overview
figure;
plot(x, z,         'b-',  'LineWidth', 1.2, 'DisplayName', 'Raw profile (flipped to peaks up)'); hold on;
plot(x, form_component, 'c--', 'LineWidth', 1.2, 'DisplayName', 'Form component (removed)');
plot(x, z_s,       'm-',  'LineWidth', 1.2, 'DisplayName', 'After form removal + S-filter');
plot(x, long_wave + mean(z), 'k--', 'LineWidth', 1.2, 'DisplayName', 'Long-wave (shown with form added back)');
plot(x, r_profile, 'r-',  'LineWidth', 1.2, 'DisplayName', 'R-profile');
xlabel('x (mm)'); ylabel('z (mm)');
title('Profile filtering overview');
legend('Location','best'); grid on; hold off;

% Plot 2: Ra bands on evaluation-length R-profile
r_mean_eval = mean(r_eval);
figure;
plot(x_eval, r_eval, 'r-', 'LineWidth', 1.5); hold on;
yline(r_mean_eval,            'k:',  'LineWidth', 1.2, 'DisplayName', 'Mean line');
yline(r_mean_eval + Ra_eval,  'r--', 'LineWidth', 1.0, 'DisplayName', '+Ra');
yline(r_mean_eval - Ra_eval,  'r--', 'LineWidth', 1.0, 'DisplayName', '-Ra');
xlabel('x (mm)'); ylabel('z (mm)');
title(sprintf('Ra bands on R-profile (Ra = %.4f mm)', Ra_eval));
legend('R-profile','Mean','+Ra','-Ra','Location','best');
grid on; hold off;

% Plot 3: Rq bands
figure;
plot(x_eval, r_eval, 'r-', 'LineWidth', 1.5); hold on;
yline(r_mean_eval,            'k:',  'LineWidth', 1.2);
yline(r_mean_eval + Rq_eval,  'g--', 'LineWidth', 1.0);
yline(r_mean_eval - Rq_eval,  'g--', 'LineWidth', 1.0);
xlabel('x (mm)'); ylabel('z (mm)');
title(sprintf('Rq bands on R-profile (Rq = %.4f mm)', Rq_eval));
legend('R-profile','Mean','+Rq','-Rq','Location','best');
grid on; hold off;

%% --- Step 17: Export processed profiles ---
exportTable = table(x, z, z_raw_centered, z_s, long_wave, r_profile, ...
    'VariableNames', {'x_mm','z_flipped_mm','z_raw_centered_mm', ...
                      'z_after_formS_mm','long_wave_mm','r_profile_mm'});
writetable(exportTable, 'ProcessedProfiles_mm.xlsx');
disp('Processed profiles exported to ProcessedProfiles_mm.xlsx');

%% --- Local functions ---
function g = isoGaussianKernel(lambda_mm, dx_mm, nPoints)
% ISO 16610-21 Gaussian profile filter, 50 percent amplitude
% transmission at cut-off lambda_c. Weighting function:
%   S(x) = (1/(alpha*lambda)) * exp(-pi*(x/(alpha*lambda))^2)
% with alpha = sqrt(ln(2)/pi) ~ 0.4697.
    alpha = sqrt(log(2) / pi);
    sigma = alpha * lambda_mm / sqrt(2*pi);

    filter_size = ceil(8 * sigma / dx_mm);
    if mod(filter_size, 2) == 0
        filter_size = filter_size + 1;
    end

    max_allowed = nPoints - 1;
    if max_allowed < 3
        error('Profile too short to build a stable Gaussian kernel.');
    end
    if filter_size > max_allowed
        filter_size = max_allowed;
        if mod(filter_size, 2) == 0
            filter_size = filter_size - 1;
        end
    end

    halfSize = (filter_size - 1) / 2;
    xg = (-halfSize:halfSize) * dx_mm;
    % ISO form scaled to sum = 1 for discrete convolution.
    g = (1/(alpha*lambda_mm)) * exp(-pi * (xg / (alpha*lambda_mm)).^2);
    g = g / sum(g);
end

function y = mirrorConv(x, kernel)
% Convolution with mirror (reflective) boundary extension.
    x = x(:); kernel = kernel(:);
    m = floor(numel(kernel) / 2);
    if m == 0
        y = x; return;
    end
    leftPad  = flipud(x(2:m+1));
    rightPad = flipud(x(end-m:end-1));
    xpad = [leftPad; x; rightPad];
    ypad = conv(xpad, kernel, 'same');
    y = ypad(m+1:m+numel(x));
end

function out = ternary(cond, a, b)
    if cond, out = a; else, out = b; end
end